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Abstract 

on : 

£NJ , We discuss hybrid Monte Carlo algorithms for odd-flavor lattice QCD sim- 

ulations. The algorithms include a polynomial approximation which enables 
us to simulate odd-flavor QCD in the framework of the hybrid Monte Carlo 
algorithm. In order to make the algorithms exact, the correction factor to the 
polynomial approximation is also included in an economical, stochastic way. 
■ We test the algorithms for n/ = 1, 1+1 and 2+1 flavors and compare results 

with other algorithms. 

o . 

1 Introduction 

-t— » . 
ctf . 

In lattice QCD simulations the standard algorithm to incorporate effects of dy- 
Qh! namical fermions is the Hybrid Monte Carlo (HMC) algorithm jl]] which is conven- 

tionally used to simulate two- flavor QCD. Due to the algorithmic limitation of the 
HMC, those simulations are limited to even numbers of degenerate flavors. In order 
• i— i , to include dynamical effects correctly, simulations of QCD with three light flavors 

(u,d,s quarks) are desirable. Simulations with an odd number of flavors can be 
performed using the R-algorithm ||] . Also for two flavors of staggered fermions [|| , 
the R-algorithm is chosen because the HMC is not applicable. The R-algorithm, 
however, is not exact: it causes systematic errors of order (At) 2 , where At is the 
step-size of the Molecular Dynamics evolution. A careful extrapolation to zero 
step-size is therefore needed to obtain exact results. Nevertheless, it is common 
practice to forego this extrapolation and to perform simulations with a single step- 
size chosen small enough that the expected systematic errors are smaller than the 
statistical ones. In this paper, we wish to emphasize that there is an alternative 
to the R-algorithm, which gives for our problem arbitrarily accurate results (for 
infinite computer time) without any extrapolation to Ar = 0@. Furthermore, we 
can make our algorithm exact with an additional, stochastic Metropolis test which 
we describe and implement. 

Liischer proposed a local algorithm, the so-called "Multiboson algorithm" ||, 
in which the inverse of the fermion matrix is approximated by a suitable Chebyshev 
polynomial. Originally he proposed it for two-flavor QCD. Borici and de Forcrand 
H noticed that the determinant of a fermion matrix can be written in a manifestly 
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positive way using a polynomial approximation, so that one can simulate odd 
flavors QCD with the multiboson method. Indeed, using this method, one-flavor 
QCD was simulated successfully The same polynomial approximation can be 
applied for the HMC |8| . The first example of a polynomial approximation within 
HMC simulations of two-flavor QCD was made by the authors of |8|, and later by 
|J. Application to odd flavors QCD is straightforward when one uses Borigi and de 
Forcrand's idea. Actually, in the development stage of Ref.pJ, one-flavor QCD was 
also simulated by HMC and it was confirmed that the two algorithmically different 
methods — multiboson and HMC — give the same distribution of plaquette values 

In preliminary work(ll]], we have simulated nj = 3 QCD by the HMC with a 
polynomial approximation, and shown that the results are consistent with those 
from the R-algorithm. However, the algorithm developed in is not exact yet. A 
polynomial approximation of moderate degree n may not be sufficiently accurate, 
especially for small quark masses. For such a case, it is important to control 
the errors coming from the polynomial approximation. In this study, in order 
to make our algorithm exact, we include the correction factor to the polynomial 
approximation into the algorithm and compare our results with those from other 
algorithms ( R-algorithm with one flavor, and two-flavor HMC ). 

A problem specific to odd flavors is the so-called " sign problem" : the fermionic 
determinant may change sign configuration by configuration and one has to include 
the sign into the statistical average (when quarks come in equal- mass pairs, the 
determinant is squared and its sign does not matter). The calculation of the sign 
of the determinant on each configuration is feasible in principle, but is very costly. 
It is common practice to assume that the QCD determinant is always positive. 
This is the case in the continuum limit, and there are indications that it is also the 
case for the lattice spacings normally considered and the quark masses currently 
reachable JT^]. Anyway, our goal here is not to study the sign problem, but rather 
to present an alternative method to the R-algorithm, which samples the absolute 
value of the determinant like the R-algorithm, but requires no extrapolation. 

Our paper is organized as follows. In Sec. 2 we describe the standard HMC 
algorithm which is used for even-flavor QCD simulations. In Sec. 3 we give 
our algorithm to simulate odd-flavor QCD, including the additional, stochastic 
Metropolis step which makes the algorithm exact. In Sec. 4-6 we show results from 
our algorithm and compare them with other algorithms. We give our conclusions 
in Sec. 7. 

2 Standard Hybrid Monte Carlo Algorithm 

The lattice QCD partition function with n/ flavors is given by 

Z = J [dU](U- f det Di)exp(-S g [U]), (1) 
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where Di is the Wilson Dirac fermion matrix, D{ = 1 — k%M\U\ with a hopping 
parameter K{ and M[C/] is the lattice Wilson hopping operator. S g [U] stands for 
the gauge action. 

The Hybrid Monte Carlo (HMC) algorithm was developed for a system con- 
taining multiples of two degenerate quark flavors. For instance, for nj = 2 flavors 
( two degenerate quarks ), one has the partition function 

Z = J[dU]detD 2 exp{-S g [U}). (2) 

Using a pseudofermion field </>, an integral representation of the determinant factor 
det D 2 can be given by 

[#t][^] exp (-0t£)t-i£)-i0) ) ( 3 ) 

where the relation, detD^ = detD, is used. The term (jy D^~ 1 D~ 1 (j) is written 
in a manifestly positive form which is essential to formulate the HMC algorithm. 
Introducing momenta Pj conjugate to the link variables Ui, the partition function 
is rewritten as 

Z = J[dU}[dP}[d^]{d<p]exp(-H), (4) 
where the Hamiltonian H is defined by 

H= l -P 2 + S g [U] + ^D^ l D- 1 ^ (5) 

where P 2 = Y.i P h 

The HMC algorithm consists of two steps: molecular dynamics (MD) evolution 
and Metropolis accept/reject step. In the MD evolution, one solves Hamilton's 
equations of motion: 

dUi dH 

m oh 

dt dUi u 

In general these equations are not solvable analytically. Usually one solves the 
equations approximately by the 2nd order leapfrog integrator which is time-reversible 
and area-preserving, thanks to which detailed balance is satisfied. After integrat- 
ing the equations, one obtains a candidate configuration (U',P') from the starting 
configuration (U, P). Next, one performs a Metropolis test with acceptance proba- 
bility p = min(l, exp(— SH)) where 5H = H{U' , P') — H(U, P). In case of rejection, 
one keeps the old U. In this study we call the HMC algorithm with Hamiltonian 
eq.(||) the standard HMC algorithm. 
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3 Algorithm for Odd Number of Flavors 

3.1 rif — 1 algorithm 

The partition function of nj = 1 QCD is given by 

Z = J[dU]detDexp(-S g [U]). (8) 

In order to simulate n/ = 1 QCD, one has to treat a single det-D which can not 
be incorporated into the standard HMC. Borici and de Forcrand || noticed that 
a single det D, when positive, can be written in a manifestly positive way when 
Liischer's polynomial approximation || is used. 

One can approximate the inverse of D by a polynomial, 

m 

l/D^P m (D) = l[(D-z k ), (9) 

k=l 

where z k are the roots of the polynomial P m (D) (a common choice is z k = 1 — 
exp(i 2iTk/{m + 1))). For a polynomial of even degree m = 2n, the roots come in 
complex conjugate pairs (z 2 k~i, z 2n + 2 - 2k ), k = Thus, eq.(||) is rewritten 

as 

n 

P 2 n{D) = l[(D- z 2k _ 1 )(D - ~z 2k _ x ) (10) 

k=l 

where z k = 1 — exp(i 2irk/(2n+l)). Using the 75 hermiticity of the fermion matrix, 
one finds that det(D — z k ) = det(D — Zky- Introducing T n (D) = ]J^ =1 (D — z 2 k-i), 
the determinant of D is written as 

det(D) = C n det(Tl(D)T n (D))-\ (11) 

with the correction factor C n given by 

C n = det(DTl(D)T n (D)) (12) 

which goes to one in the limit n — > 00. An integral representation of det(T^(D)T n (D)) 
is given by 

det(T ? t( J D)r n ( J D))- 1 ~ J [#t][ # ] exp(-0t^(Z))r n ( J D)0). (13) 

Note that the term 4>^T^(D)T n (D)(f) is Hermitian positive, which can not be real- 
ized for one flavor in the standard HMC formulation. 

Now, the nj = 1 Hamiltonian for our HMC algorithm can be defined as 

H n f =l = lp2 + Sg[u] + 0t T t (jD)Tn(jD ) . ( 14 ) 
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Making use of this Hamiltonian one can perform HMC simulations for rif = 1 
QCD. 

The Hamiltonian defined by eq.(|l4|) is an approximation to the exact one, which 
generates configurations with the measure ~ det(T T [(D)T n (D)) _1 . The quality of 
this approximation is measured by the deviation of C n (eq.(|i~2|)) from one. To 
make the algorithm exact, one has to include the correction factor C n into the 
algorithm. There are several possibilities to incorporate C n in the measure fl3|| . 
One possibility is to make a global Metropolis test with probability 



P 



mm 



Cn 1 

C n 



(15) 



where C' n = C n [U'] and C n = C n [U] correspond to a new candidate configuration 
V and a starting configuration U respectively. Since the correction factor is the 
determinant of a large matrix, its direct calculation is not feasible 1 . An econom- 
ical way is to form an unbiased estimator of the ratio C' n /C n and to use a noisy 
Metropolis test M, PTJ]. The same correction factor appears in the nj = 1 multibo- 
son algorithm |Q]7where the ratio was estimated by rewriting the correction factor 
using another high-quality polynomial T r (D) (r » rt)|7], 16] 



as 



Cn 



Hm detT n (D)T n (D) 
r ^°° detT ] r {D)T r (D) 



(16) 



Using this form of the correction factor, the ratio C' n /C n is written as 



C'n 

C-n. 



det 



X*X' 



X^X 
J drfdrj exp( 



iX'^Xrif 



J drfdr] exp(— \r]\ 2 ) 
<exp(-|X'- 1 X7 ? | 2 + |r ? | 2 ) > v 



(17) 



where X = T n (D)T r (D)^ 1 . Thus the ratio can be estimated by calculating 
exp(— \X'~ 1 Xr]\ 2 + \rj\ 2 ) with only one Gaussian random vector r\. However con- 
vergence is slow, and there may be a technical difficulty using a high-quality poly- 
nomial: the above estimation includes a number of " matrix(D) x vector + vector 
" type calculations, which results in divergence for high-degree polynomials due to 
the roundoff errors of computers [jO]]. 

We solve this problem very economically here, by proposing to estimate the 
ratio C' n /C n from unbiased estimators of {C' n /C n ) 2 . 

(C n /C n ) 2 is easily estimated by 



C' n 

C n 



det 



D' 2 P 2n {D^P 2n {D') 
D 2 P 2n (DyP 2n (D) 



det 



1 In jjij the correction factor was computed exactly by the Lanczos method. This approach is 
limited to small lattices. 
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J drfdr] exp(— \t]\ 2 ) 
= < expi-lW'^Wvl 2 + \rj\ 2 ) > v (18) 

where W = DP2 n (D). Then the problem is to estimate (C' n /C n ) using unbiased 
estimators of (C' n /C n ) 2 only. This can be accomplished by evaluating stochasti- 
cally the Taylor expansion of y/x as shown in WR for e x . Expanding about x = 1, 
one writes 

= 1 + \ e + E CfeR^ (19) 

fc=2 

where = 2k 2 ~)/' c k-i £]0, 1[. First, the left-hand side is assigned value 1. Next, 
with probability 5 one computes a stochastic estimator of ei = x — 1 and adds 
it to the left-hand side. Then, with probability 2 ^^T 3 Cfc— 1 1 fc=2 = \ one computes 
a second stochastic estimator of £2 and adds — e\€2 to the left-hand side, and so 
on. In our case, x = {C' n /C n ) 2 . When the procedure terminates, one obtains an 
unbiased estimator of yfx. There is a difficulty if the estimator becomes negative: 
it cannot be used as a probability in the Metropolis test (eq.([l5|)). Such positivity 
violations can be reduced by increasing the degree n, to the point where they are 
never observed during the whole simulation. 
Our algorithm is summarized as follows. 

1. HMC: we perform Molecular Dynamics with the approximate Hamiltonian 
eq.(|i~4|) and obtain a candidate configuration U' . Then we do a Metropolis 
test with acceptance probability p = min(l, 

exp (-£#)) This 

removes the 

step-size integration error. 

2. If the candidate configuration is accepted, we estimate the ratio C' n /C n 
stochastically by using the method explained above. 

3. We perform another Metropolis test with acceptance probability p = min(l, y/x), 
where y/x is an unbiased estimator of C' n /C n . This removes the error of the 
polynomial approximation to detD. 

4. If accepted, we take U' as a new configuration. Otherwise we keep the old 
configuration. 

This algorithm samples the measure Vdet 2 D = \detD\, just like the R- 
algorithm, but does so in a more efficient way which avoids the extrapolation 
to zero step-size of the latter. For very small quark masses in the Wilson formu- 
lation, the Dirac determinant may become negative for some background gauge 
fields. This happens for a small fraction of the gauge ensemble, which goes to 
zero in the continuum limit. Therefore, sampling the measure |detD| does not 

2 Here, as well as in the HMC algorithm, one could use the Glauber acceptance probability 
p = (1 + cxp(— SH))' 1 instead of the Metropolis one. This does not seem to be a judicious choice 
however, since even for an exact MD evolution (8H = 0) the acceptance will only be 1/2. 
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affect the continuum limit of the lattice results. Nevertheless, it is possible, if de- 
sired, to correct for negative determinants by multiplying the contribution of the 
corresponding configurations by — 1 in the Monte Carlo ensemble. This requires 
identifying possible negative real Dirac eigenvalues in the configurations of that 
ensemble, which can be done e.g. with the Arnoldi method. 

Note that the choice of approximating polynomial is in principle arbitrary: 
the polynomial approximation error is removed at Step 3 above. However, a poor 
choice of polynomial, e.g. whose domain of approximation does not cover the com- 
plete spectrum of the Dirac operator, will result in difficulty maintaining positivity 
of the estimator of C' n /C n , and in longer autocorrelation times. 

3.2 rif = 2 + 1 algorithm 

The partition function of nf = 2 + 1 QCD is given by 

Z = J [dU] det D\ det D t exp(-S g [U]). (20) 

This system consists of two degenerate quark flavors with a hopping parameter n k 
and one flavor with k\ . An integral representation of the determinant factor det Di 
is obtained using the polynomial approximation as in eq. (|l3j) , whereas det-Dj? is 
expressed as in eq.(||) 

det Di ~ det- 1 (Tt(A)T n (A)) ~ J exp(-^(A)r„(A)0O. ( 21 ) 

and 

det Dl ~ J MWk] expC-^^I -1 ^ Vfc). (22) 
Combining eq.(pTI) and (22) one can define the following Hamiltonian, 



H n f= 2+1 = l -P 2 + S g [U] + <\>\D\- x D k -H k + ^Tt(A)T n (A¥z. (23) 
Alternatively, one can express the determinant factors using only one vector <f> as 

det^detA- y[d < /» t ][#]exp(- ( /.t^r 1 r^(A)r n (A) J Dfc V)- (24) 

Then one obtains another Hamiltonian 

H n f= 2+1 = \P 2 + S g [U] + ^Dl'^iD^iD^cp. (25) 

Both definitions of Hamiltonian can be used for HMC simulations. In this study 
we take eq. (|23|) . 

The correction factor to the rif = 2 + 1 Hamiltonian is given by 

C n = det(ATt(A)T n (A))- (26) 
This factor can be included in the same way as explained for the rif = 1 algorithm. 
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3.3 rif = 1 + 1 algorithm 

The partition function of n/ = 1 + 1 QCD is given by 

Z = J [dU] det D k det D l exp(-S g [U]). (27) 

This system consists of two non-degenerate quark flavors with hopping parameters 
kj~ and Ki- Each determinant factor can be expressed in terms of pseudofermion 
fields using polynomial approximations as 

det D k ~ Jid^d^expi-^lT^D^TniDk)^), (28) 

detA~ J[d<j)i][d4}eM-4 T l( D i) T rn{Di)<t>i). (29) 
Then one defines the rif = 1 + 1 Hamiltonian as 

# n/=1+1 = ^p 2 + 5 3 [?7] + <j>lTl{D k )T n {p k )^ h + 0fr^(A)r m (A)^ (30) 

Using one field only, we can express the two determinant factors at once as 



det£> fc detA~ J [#][# t ]exp(-^Tt( J D fc )Tt 4 (A)T m (A)r n ( J D fc )</>). (31) 
This expression results in the following Hamiltonian, 

H n f= 1+1 = X -P 2 + S g [U] + t^T^(Dk) r r^ l (Di)T m (Di)T n (Dk)<l>. (32) 

In this study we use the definition of eq.([30|) for HMC simulations. The correction 
factor of nt = 1 + 1 Hamiltonian is given by 



Cnm — det 



D k Tl(D k )T n (D k ) ■ A7^(A)T m (A)l • (33) 



4 Numerical results for n/ = 1 + 1 

We have tested the above algorithms, using relatively small lattices in order to 
obtain results of sufficient accuracy to expose tiny systematic errors caused by a 
low-order polynomial approximation or a finite step-size (for the R-algorithm) . 

We first show results of rif = 1 + 1 QCD with degenerate quark masses. This 
allows a direct comparison of our results with those of standard rif = 2 HMC. The 



algorithm we use here is the one discussed in Sec. 3.3, where we set We 
use an 8 lattice at = 5.30 and k = 0.156,0.158 and 0.160. Boundary conditions 
are anti-periodic in all directions. Since the two quarks have the same mass, we use 
identical polynomials with the same degree n for the polynomial approximation, 
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ie. T n (D) = T m (D) and n = m and write the correction factor as C n = C nm . We 
measure lxl, 1x2 and 2x2 Wilson loops and vary the degree n of T n (D). The 
average values are taken over about 2000-8000 trajectories. Details are shown in 
Table ||. Average values of Wilson loops are displayed in Figs.|^-^ with those from 
the HMC. Results of Wilson loops are consistent with those from the standard 
HMC except for small n. These results suggest that we do not need to take a very 
large n. Quantitatively, however, from results of Wilson loops only it is not clear 
which n should be chosen. 

Fig.[i] shows the average value of \C' n /C n — 1| as a function of the degree n. 
Negative values of C' n /C n observed for small n are not included in this average. 
Except for small n, C' n /C n converges to one exponentially as n increases. The 
plateau seen at small n is due to the following reason. When the polynomial 
approximation is inaccurate, i.e. for small n, |VF /_1 P^r/| 2 takes a large value, and 
the estimated value of (C' n /C n ) 2 via eq.(|l^) will always be very small. If we apply 
eq.(||) with {C' n /C n ) 2 = (i.e. e = -1), we obtain (C' n /C n ) = 0.429 when the 
average excludes negative values. This value is consistent with our results. Anyhow 
we are not interested in such a small n and should increase n until no negative 
value of C' n /C n is observed. 

Fig.|5] shows the positivity- violation (PV) rate of C' n /C n . The PV appears to 
be suppressed exponentially as n increases and one can easily choose the degree n 
such that no PV would be observed within desired statistics. Within our statistics, 
no negative value of C' n /C n was observed for n > 22, 24 and 30 at k = 0.156, 0.158 
and 0.160 respectively. However, these values of n will depend on statistics: with 
high-statistics one may observe a small number of negative values even at large n. 
In order to remove this dependency, let us fix the rate of positivity violation. In 
this study we are dealing with statistics of O(10 3 — 10 4 ) trajectories. So we set the 
PV level to 10 -4 . From Fig.0, roughly speaking, we find that this level corresponds 
to n 23,28,34 for k = 0.156, 0.158 and 0.160 respectively. The lines indicated 
by "No Positivity Violation" in Figs.[l|-|| correspond to these values of n. In turn, 
as seen in Fig|| these values of n correspond to an acceptance ~ 95% in Step 3 
of our algorithm. This suggests that typically one needs an acceptance of 95% or 
higher to maintain positivity of C' n /C n at the level of 10 -4 . In terms of C' n /C n , 
roughly speaking, 95% acceptance corresponds to < \C' n /C n — 1| 0.1. Under 
these conditions, the contribution of configurations with negative C' n /C n to the 
total ensemble, if present at all, will be well below the statistical fluctuations, and 
the way they are dealt with is unimportant. We simply reject such configurations. 



5 Numerical results for nj = 1 

We take a 6 4 lattice at (3 = 5.45 and n = 0.160 with periodic boundary condi- 



tions in all directions. Simulations are performed using the algorithm of Sec. 3.1 
Average values of Wilson loops are shown in Fig.|?] together with results from the 
R-algorithm extrapolated to zero step-size. The average values are taken over 
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10000-36000 trajectories. See Table |3] for details. Our results are consistent with 
those from the R-algorithm except for small n. 

Fig|| shows the average value of \C' n /C n — 1|. Again, we see exponential con- 
vergence as n increases, except for the plateau at small n. 

Fig.|9] shows the PV rate of C' n /C n . Positivity violation appears to be sup- 
pressed exponentially as n increases. Within our statistics, positivity was main- 
tained for n > 30. As discussed in Sec.|4|, adopting a PV level of 10~ 4 , we find that 
this level corresponds to n 33. 

Fig.|l^ shows the average acceptance of the Metropolis test Step 3. The PV level 
of 10~ 4 corresponds again to about 95% acceptance, and to < \C' n /C n — 1| 0.07. 



6 Numerical results for rif = 2 + 1 

We take a 6 lattice at f3 = 5.30 and k = 0.156 with periodic boundary conditions 
in all directions. Simulations are performed using the algorithm of Sec. [3.2| . Setting 
k = k/j = ki, we can simulate n/ = 3 QCD and compare results with those from 
rif = 3 R-algorithm. The average values are taken over 9000-16000 trajectories. 
See Table || for details. Our results are consistent with those from the R-algorithm 
extrapolated to zero step-size, except for small n. 



Fig. 12 shows the average value of \C' n /C n — 1|. Again we observe exponential 



convergence as n increases, except for small n. 



Fig. 13 shows the positivity-violation (PV) rate of C' n /C n . Within our statistics, 
positivity was maintained for n > 24. As discussed in SecQ, adopting a PV level 
of 10 -4 , we find that this level corresponds to n ~ 27. 



Fig.|14j shows the average acceptance of the Metropolis test Step 3. A PV level 

-4 



of 10 corresponds to about 95% acceptance and to < \C' n /C n — 1| 0.1. 



7 Conclusions 

We have discussed HMC algorithms for odd-flavor QCD simulations, based on a 
polynomial approximation of the inverse Dirac operator D~ l . The algorithms can 
be made exact with a correction factor which can be easily incorporated in an 
additional, stochastic Metropolis test. We have tested the algorithms for nj = 
1,1 + 1 and 2+1 flavors. The results are consistent with those from the standard 
HMC and R-algorithm. 

The estimator of the correction factor C' n /C n should be positive if it is to be 
used in a Metropolis test. We observe that positivity violation is suppressed expo- 
nentially as n increases. Therefore, one can choose in advance a value of n which 
will suffice to avoid negative correction factors in the simulation for the desired 
statistics. When one fixes the positivity-violation level to 10 -4 , this corresponds 
to about < \C' n /C n — 1| >ss 0.07 — 0.1. Equivalently, this level corresponds to 
« 95% acceptance in the Metropolis test with p = min(l, C' n /C n ). 
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Table 1: Step-size At, number of trajectories and acceptances for n/ = 1 + 1 QCD 
simulations on an 8 4 lattice at (3 = 5.30. Acc.(HMC) stands for the acceptance 
with p = min(l, exp(— 5H)) which corrects the integration step-size error, and 
Acc.(CF) is the acceptance from the correction factor with p = min(l, C' n /C n ), 
which corrects the error of the polynomial approximation. 



K 


n 


At 


trajectories 


Acc.(HMC)% 


Acc.(CF)% 


0.156 


4 


0.05 


3000 


78 


52 




fi 

\j 


n D5 


3nnn 


sn 

OVJ 


51 




8 


0.05 


3000 


80 


51 




10 


0.05 


4500 


80 


49 




20 


0.05 


3000 


81 


88 




22 


0.05 


4000 


81 


93 




24 


0.05 


3000 


80 


96 




30 


0.05 


4000 


81 


98.9 




40 


0.05 


4200 


82 


100 




46 


0.05 


3600 


80 


100 


0.158 


4 


0.05 


3000 


78 


49 




6 


0.05 


4500 


77 


50 




8 


0.05 


3000 


80 


50 




10 


0.05 


4500 


80 


49 




20 


0.05 


5000 


81 


77 




24 


0.05 


4000 


80 


90 




30 


0.05 


4000 


80 


98 




40 


0.05 


4900 


80 


99.7 




46 


0.05 


4500 


80 


99.9 


0.160 


4 


0.05 


5000 


78 


51 




6 


0.05 


4000 


79 


50 




8 


0.05 


7500 


80 


50 




10 


0.05 


4500 


80 


49 




20 


0.05 


6000 


79 


61 




24 


0.05 


3600 


80 


77 




30 


0.05 


3200 


81 


92 




40 


0.05 


2800 


81 


98.8 




46 


0.05 


3000 


81 


99.7 
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Table 2: n/ = 1 + 1 simulation results for lxl, 1x2 and 2x2 Wilson loops on 



an 8 la 


,tice at 


P = 5.30. 






K 


n 


1 x 1 


1x2 


2x2 


0.156 


4 


0.48846(37) 


0.25172(48) 


0.07455(37) 




6 


0.48637(28) 


0.24927(34) 


0.07276(22) 




8 


0.48502(44) 


0.24780(51) 


0.07176(39) 




10 


0.48566(48) 


0.24845(60) 


0.07235(45) 




20 


0.48548(39) 


0.24824(45) 


0.07212(32) 




22 


0.48568(16) 


0.24853(20) 


0.07234(16) 




24 


0.48589(34) 


0.24876(42) 


0.07246(32) 




30 


48566(251 


24845 ( 311 


07224(231 




40 


0.48587(31) 


0.24877(40) 


0.07245(31) 




46 


0.48613(28) 


0.24906(36) 


0.07271(26) 




HMC 


0.48594(33) 


0.24883(41) 


0.07256(32) 


0.158 


4 


0.49274(48) 


0.25698(58) 


0.07843(46) 




6 


0.49129(44) 


0.25531(57) 


0.07730(48) 




8 


0.49032(48) 


0.25403(55) 


0.07628(37) 




10 


0.49055(46) 


0.25431(57) 


0.07646(40) 




20 


0.49020(24) 


0.25391(28) 


0.07615(19) 




24 


0.49055(29) 


0.25439(35) 


0.07661(26) 




30 


0.49003(39) 


0.25373(47) 


0.07603(36) 




40 


0.49027(42) 


0.25401(52) 


0.07634(40) 




46 


0.49009(21) 


0.25378(25) 


0.07608(20) 




HMC 


0.49026(20) 


0.25397(25) 


0.07624(19) 


0.160 


4 


0.49814(42) 


0.26347(52) 


0.08326(42) 




6 


0.49604(48) 


0.26110(61) 


0.08180(49) 




8 


0.49546(55) 


0.26030(69) 


0.08102(54) 




10 


0.49515(34) 


0.25982(45) 


0.08054(40) 




20 


0.49457(34) 


0.25923(42) 


0.08016(34) 




24 


0.49569(29) 


0.26055(35) 


0.08120(30) 




30 


0.49555(33) 


0.26038(49) 


0.08102(34) 




40 


0.49564(41) 


0.26051(51) 


0.08106(40) 




46 


0.49529(43) 


0.26014(55) 


0.08089(42) 




HMC 


0.49573(48) 


0.26068(62) 


0.08131(49) 
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Table 3: Same as in table [T], but for rif = 1 QCD simulations on a 6 4 lattice at 
P = 5.45. 



n 


At 


trajectories 


Acc.(HMC)% 


Acc.(CF)% 


2 


0.04 


5000 


93 


49 


4 


0.04 


22000 


93 


50 


6 


0.04 


36000 


92 


50 


8 


0.04 


20000 


93 


49 


10 


0.05 


15000 


89 


50 


20 


0.04 


14000 


93 


75 


24 


0.04 


10500 


93 


86 


30 


0.04 


14300 


93 


94 


34 


0.05 


16500 


89 


97 


40 


0.05 


15000 


89 


99 


46 


0.05 


22000 


89 


99.4 



Table 4: rif = 1 simulation results for lxl, 1x2 and 2x2 Wilson loops on a 6 4 
lattice at (3 = 5.45. Results from the R-algorithm were extrapolated to step-size 
At = 0, based on simulation results at At = 0.0125, 0.025, 0.0333333 and 0.04. 



K n 


1 x 1 


1 x 2 


2x2 


0.160 2 


0.49637(60) 


0.25899(73) 


0.07794(59) 


4 


0.51235(54) 


0.27934(70) 


0.09420(60) 


6 


0.51305(51) 


0.28041(67) 


0.09543(57) 


8 


0.51199(36) 


0.27896(46) 


0.09412(40) 


10 


0.51206(36) 


0.27904(47) 


0.09422(40) 


20 


0.51195(47) 


0.27885(61) 


0.09399(53) 


24 


0.51174(39) 


0.27863(51) 


0.09387(45) 


30 


0.51206(41) 


0.27906(53) 


0.09421(47) 


34 


0.51207(31) 


0.27907(40) 


0.09423(36) 


40 


0.51190(28) 


0.27886(35) 


0.09402(30) 


46 


0.51143(33) 


0.27822(43) 


0.09347(37) 


R-algorithm 


0.51163(29) 


0.27850(38) 


0.09377(35) 
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Table 5: Same as in table |T| but for nj = 2 + 1 QCD simulations on a 6 4 lattice at 
P = 5.30. 



n 


At 


trajectories 


Acc.(HMC)% 


Acc.(CF)% 


2 


1/18 


9000 


80 


50 


4 


1/18 


12100 


80 


51 


6 


1/18 


15700 


81 


51 


8 


1/18 


12000 


80 


50 


10 


1/18 


9000 


79 


51 


20 


1/18 


12500 


81 


82 


24 


1/18 


10000 


81 


92 


30 


1/18 


12500 


81 


97 


40 


1/18 


14000 


81 


99.5 


46 


1/18 


16000 


80 


99.9 



Table 6: rif = 2 + 1 simulation results for lxl, 1x2 and 2x2 Wilson loops 
on a 6 4 lattice at f3 = 5.30. Results from the R-algorithm were extrapolated to 
step-size At = 0, based on simulation results at At = 0.0125, 0.0333333, 0.04 and 



0.05. 



k n 


1 x 1 


1 x 2 


2x2 


0.156 2 


0.50032(46) 


0.26559(58) 


0.08444(48) 


4 


0.51785(77) 


0.2882(10) 


0.10343(95) 


6 


0.5206(11) 


0.2919(15) 


0.1071(14) 


8 


0.52072(92) 


0.2922(12) 


0.1075(11) 


10 


0.52005(69) 


0.29124(92) 


0.10648(85) 


20 


0.52161(71) 


0.29338(97) 


0.10859(94) 


24 


0.51993(77) 


0.2911(10) 


0.10639(99) 


30 


0.52131(98) 


0.2930(13) 


0.1082(12) 


40 


0.52009(75) 


0.2913(10) 


0.10657(98) 


46 


0.51950(79) 


0.2905(11) 


0.1059(10) 


R-algorithm 


0.5204(10) 


0.2917(14) 


0.1069(14) 
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Figure 1: rif = 1 + 1 results of Wilson loop on an 8 4 lattice at (3 = 5.30 and 
k = 0.156. 
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Figure 2: Same as in fig.| but for k = 0.158. 



19 



8 P=5.30 k=0.160 



0.5 
0.499 
0.498 
: 0.497 
: 0.496 
0.495 
0.494 
0.493 







10 



o n f = 1+1 with CF 
□ HMC 



No Positivity Violation 



20 



30 

n 



40 



50 



60 



0.265 
0.264 
0.263 



cd 0.262 



: 0.261 
0.26 
0.259 
0.258 







10 



8 P=5.30 k=0.160 



20 



o n f = 1+1 with CF 
□ HMC 



No Positivity Violation 



30 

n 



40 



50 



60 



8 0=5.30 k=0.160 



0.084 



0.083 



A 

<£! 0.082 

CM 



0.081 



0.08 



0.079 







10 



o n f = 1+1 with CF 
□ HMC 



No Positivity Violation 



20 30 40 50 

n 



60 



Figure 3: Same as in fig.@ but for k = 0.160. 
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Figure 4: \C' n /C n — 1| as a function of degree n for nj = 1 + 1 at k = 0.156, 0.158 
and 0.160. R in the figure stands for C' n /C n . C' n /C n are estimated stochastically 



by eq.(18) and negative values are not used in the average. 
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Figure 5: Rate of positivity violation of C' n /C n for rif 
degree n. 
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Figure 6: Acceptance of the Metropolis test with p = min(l, C' n /C n ) for nj = 1 + 1. 



22 



6 4 0=5.45 k=0.160 



0.514 



0.513 



V 0.512 



0.511 



10 



o n f =1 with CF 
□ R-algorithm 



No Positivity Violation 



20 



30 

n 



40 



50 



60 



0.282 

0.281 

A 0.28 
x 

V 0.279 
0.278 
0.277 



10 



6 |3=5.45 k=0.160 



o n f =1 with CF 
□ R-algorithm 



No Positivity Violation 



20 



30 

n 



40 



50 



60 



0.097 



0.096 



A 

$ 0.095 



0.094 



0.093 



0.092 



6 0=5.45 k=0.160 



o n f =1 with CF 
□ R-algorithm 



No Positivity Violation 



10 20 30 40 50 60 



Figure 7: rif = 1 results of Wilson loop on a 6 4 lattice at /3 = 5.45 and k = 0.160. 
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Figure 9: Rate of positivity violation of C' n /C n for nj 
n. 
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Figure 10: Acceptance of Metropolis test with p = min(l, C' n /C n ) for rif = 1. 
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Figure 11: nj = 2 + 1 results of Wilson loop on a 6 4 lattice at (3 = 5.30 and 
k = 0.156. 
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Figure 12: | C' n /C n — 1| for nj = 2 + 1 as a function of degree n. R in the figure 
stands for C' n /C n . 
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Figure 13: Rate of positivity violation of C' n /C n for rif = 2 + 1 as a function of 
degree n. 
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Figure 14: Acceptance of Metropolis test with p = min(l, C' n /C n ) for rif = 2 + 1. 
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